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ABSTRACT 


An analysis and comparison of daily and yearly solar irradiation from the satellite CM SAF database and a 
set of 301 stations from the Spanish SIAR network is performed using data of 2010 and 2011. This analysis 
is completed with the comparison of the estimations of effective irradiation incident on three different 
tilted planes (fixed, two axis tracking, and north-south horizontal axis) using irradiation from these two 
data sources. Finally, a new map of yearly values of irradiation both on the horizontal plane and on inclined 
planes is produced mixing both sources with geostatistical techniques (kriging with external drift, KED). 

The Mean Absolute Difference (MAD) between CM SAF and SIAR is approximately 4% for the irradiation 
on the horizontal plane and is comprised between 5% and 6% for the irradiation incident on the inclined 
planes. The MAD between KED and SIAR, and KED and CM SAF is approximately 3% for the irradiation on 
the horizontal plane and is comprised between 3% and 4% for the irradiation incident on the inclined 
planes. 


CM SAF The methods have been implemented using free software, available as supplementary material, and the 

SIAR data sources are freely available without restrictions. 
© 2013 Elsevier Ltd. All rights reserved. 
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Nomenclature 


AEMET Spanish Meteorology Agency 

By estimated coefficients of the deterministic model in 
kriging with external drift 

BSRN Baseline Surface Radiation Network 

CM SAF Satellite Application Facility on Climate Monitoring 


AG(0) difference between G™^ (0) and G>“*(0) 
AGWSAF(Q) difference between G? (0) and G™S4F(Q) 


AGpan(0) difference between G'®” (0) and G" (0) 


G™SAF(Q) yearly global irradiation on the horizontal plane data 


obtained by estimations from CM SAF 
yearly global irradiation on the horizontal plane data 
obtained by on-ground measurements from SIAR 


AGH RED difference between G? (0) and G™ (0) 


AGH KED difference between G"(0) and G” (0) 


GAR (0) 


AGer difference between G^" and Gé"* 
€(Sg) interpolated residual in kriging with external drift 


GEP(0) yearly global irradiation on the horizontal plane 


estimated with kriging with external drift 


Gy (0) yearly effective global irradiation on the inclined 
plane estimated with kriging with external drift 

y(h) semivariogram function 

>ch) estimator of the semivariogram function 


Govs4F yearly effective global irradiation incident on differ- 
ent planes estimated from data from CM SAF 

GJ yearly effective global irradiation incident on differ- 
ent planes estimated from data from SIAR 

h separation vector between two locations 

KED kriging with external drift 

Ài kriging weights determined by the spatial depen- 


dence structure of the residual 
LUT look-up table 


(So) fitted deterministic part of the random spatial field at 
a new location 

MAB Mean Absolute Bias 

MAD Mean Absolute Difference 

MBD Mean Bias Difference 


OK ordinary kriging 


qk(Sọ) auxiliary predictors obtained from the fitted values of 
the explanatory variable at the new location in kri- 
ging with external drift 

RMSD _ biased Root Mean Square Difference 

RMSD* unbiased Root Mean Square Difference 

RTM radiative transfer model 

SIAR Agroclimatic Information System for Irrigation 

SIS shortwave incoming solar radiation 

Zz kriging estimation of the random spatial field 

Z(S) random spatial field 


1. Introduction 


Nowadays, with a wide range of applications in agriculture, 
climate monitoring and renewable energies, research in solar 
irradiation is a very demanded field. 

Solar irradiation can be evaluated by processing images from 
satellites or by on-ground measurements with pyranometers in 
meteorological stations. The high cost of these meteorological 
stations and the requirement of specific and periodic calibrations 
explain the low density of the existing networks in many countries, 
although this kind of measurements is reliable to elaborate solar 
irradiation maps [1]. The satellite models need to be validated and 
refined with high quality measurements, which are provided by on- 
ground stations [2]. The satellite estimates present a wide spatial 
and temporal coverage, but their spatial resolution is in the range of 
kilometres, which in many applications may not be sufficient and 
this can be improved with geostatistics [3]. The high degree of site 
dependence of solar irradiation makes geostatistics suitable to 
evaluate the spatial distribution of solar irradiation and to build 
maps with pyranometers measurements [4,5]. 

Geostatistics were firstly applied in the study and estimation 
of ore resources [6], soil properties [7] and afterwards, in fields 
such as on-ground water analysis [8] and solar irradiation maps 
with kriging techniques [3,9]. Residual and ordinary kriging have 
been applied to elaborate solar irradiation maps taking into 
account elevation and cloudiness as significant variables [10], or 
topographic shadow cast and elevation [11], and also with 
artificial neural networks (ANN) with temperature and precipita- 
tion as inputs [12]. Kriging with external drift (KED) has been 
useful to develop solar irradiation maps using multiple linear 
regression (MLR) models [13]. Comparing solar irradiation maps 
obtained with different techniques and inputs is necessary to 
assess the divergence of the estimates. The MESOR project 
compared EnMetSol, Helioclim-2, NASA SSE version 6, Satel- 
Light and SOLEMI databases obtained with satellite estimates 
and ESRA, PV GIS Europe, and Meteonorm version 6.1 databases 


generated from geostatistical models and meteorological obser- 
vations in Europe [1]. 

Recently, the Spanish Agency of Meteorology (AEMET) has 
released a new solar irradiation atlas for Spain (the former 
was of 1984) [14] providing monthly, seasonal and annual 
average of global, direct and diffuse irradiation on the hor- 
izontal plane with a resolution of 3 km using monthly data 
sets from 1983 to 2005 of CM SAF. Besides, a validation 
process has been developed comparing CM SAF data with 
uninterrupted registers from 2003 to 2005 of 29 meteorolo- 
gical stations from the National Radiometric Network (RRN) of 
AEMET. On the other hand, for direct irradiation, only two 
ground stations, with uninterrupted data from 1992 to 2005, 
were selected. The Mean Absolute Deviation (MAD) obtained 
from this validation process for global monthly average is 
12.23 W/m? (6.7%), which is slightly higher than the CM SAF 
target of 10 W/m7. It is important to underline that the AEMET 
global irradiation atlas is restricted to the horizontal plane. 

This paper innovates with an analysis and comparison of solar 
irradiation from the CM SAF database (Section 2.1) and a large set 
of stations, considering 301 meteorological stations (versus the 29 
of the aforementioned assessment by AEMET) from the Spanish 
SIAR network (Section 2.2), and with the estimation of effective 
irradiation incident on three different tracking planes. 

Therefore, the contribution of this paper is threefold: 


e Analysis and comparison of daily and yearly global irradiation 
on the horizontal plane obtained by on-ground measurements 
and satellite estimate data. 

e Analysis and comparison of yearly global irradiation incident 
on different tilted planes (fixed, two axis tracking on azimuth 
and solar elevation, north-south horizontal axis) estimated 
from these two data sources. 

e Elaboration of a new map of yearly values of irradiation both 
on the horizontal plane and on inclined planes with a smooth 
combination of both sources using geostatistical techniques. 
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Fig. 1. Hovmöller plot with the time evolution of the daily horizontal irradiation 
(Wh/m?) as published by CM SAF, averaged along 10°W to 5°E. 


The analysis comprises daily irradiation data of 2010 and 
2011. The global irradiation on the horizontal plane is compared 
both in a daily and a yearly basis, while the effective irradiation 
incident on different planes is only examined in a yearly basis. In 
order to ease the discussion of results, the yearly analysis is 
carried out with the averages of 2010 and 2011. 

To enable reproducible research [15], the methods have been 
implemented using free software (Section 6). Both the source code 
and the data sources are freely available without restrictions. 


2. Radiation data sources 
2.1. CM SAF 


The Satellite Application Facility on Climate Monitoring (CM 
SAF) [16] is a joint venture of the Royal Netherlands Meteorolo- 
gical Institute, the Swedish Meteorological and Hydrological 
Institute, the Royal Meteorological Institute of Belgium, the 
Finnish Meteorological Institute, the Deutscher Wetterdienst, 
Meteoswiss, the UK MetOffice, with the collaboration of the 
European Organization for the Exploitation of Meteorological 
Satellites (EUMETSAT). The CM SAF was funded in 1992 to 
retrieve, archive, and distribute climate data to be used for 
climate monitoring and climate analysis. The spatial resolution 
of the different products ranges from 15 km? to 90 km? [17]. 

The CM SAF provides two categories of data: operational 
products and climate data. The operational products are built on 
data that is validated with on-ground stations and then is 
provided in near real time to develop variability studies in diurnal 
and seasonal time scales. However, climate data are long-term 
data series to assess inter-annual variability [18]. 

In this study, the shortwave incoming solar radiation product 
(SIS) is selected with a spatial resolution of 15 km?, available as 
daily and monthly averages (Figs. 1 and 2). SIS collates shortwave 
radiation (0.2—4 um wavelength range) reaching an horizontal 
unit earth surface obtained by processing information from 
geostationary satellites (SEVIRI sensor on board of the METEOSAT 
Second Generation (MSG)) and also from polar satellites (AVHRR 
sensor on NOAA polar satellites) [17] and then validated with 
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Fig. 2. Average of yearly horizontal irradiation (kWh/m?) on the horizontal plane 
as published by CM SAF during 2010 and 2011. 
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Fig. 3. Meteorological stations of the SIAR network. The color key indicates the 
altitude (m). Those stations whose average yearly absolute difference from the CM 
SAF values is higher than 5% are displayed with red points. (For interpretation of 
the references to color in this figure caption, the reader is referred to the web 
version of this article.) 


high-quality on-ground measurements from the Baseline Surface 
Radiation Network (BSRN).! 

In this paper, SEVIRI data has been selected following the CM 
SAF recommendation of these data to be used for latitudes south- 
ern 65°N [19]. Validation of SEVIRI SIS data with 4 BSRN stations 
showed that more than 90% of the values are below the accuracy 
target value of 10 W/m? (plus the uncertainty of the ground based 
measurements). Besides, the absence of a trend in the bias 
demonstrates the stability and homogeneity of the product [20]. 

The method for retrieving the solar surface irradiance 
employed by CM SAF is based on the libRadtran radiative transfer 
model (RTM) [22] in combination with a new approach of several 
parameterizations and eigenvector look-up tables (LUT). A LUT is 
a data structure with discrete pre-computed RTM results for a 
variety of atmospheric and surface states. Thus, the surface 
irradiance (transmittance multiplied by extraterrestrial incoming 
solar flux density) for a given atmospheric state can be obtained 
by interpolation, through the LUTs, for each satellite pixel and 
time. Therefore, with a LUT approach the results are similar to 
those obtained with a RTM reducing computation costs [23]. 


1 http://www.bsrn.awi.de/en/home/. 


F. Antonanzas-Torres et al. / Renewable and Sustainable Energy Reviews 21 (2013) 248-261 


GEMSAF (0) 


á Coan) 


kriging 
with 

external 

drift 


T 
fm 


251 


C> 


Effective 


Fig. 4. Organization of the analysis procedure. Ellipses represent point data sets (values from the meteorological stations, for example) and rectangles denote raster maps 
(values from CM SAF, for example). The red color is used to identify the original sources, green for comparison results, and blue for transformation results (geostatistical 
interpolation or effective irradiation). (a) Horizontal irradiation comparison. (b) Effective irradiation comparison. (For interpretation of the references to color in this figure 


caption, the reader is referred to the web version of this article.) 


The CM SAF method still can be improved by a better semi- 
empirical adjustment of cloud effects and by improved meteor- 
ological information about aerosols and snow cover maps [23]. 
In fact, one main goal of the Continuous Development and Operations 
Phase of the CM SAF (2007-2012) is to improve all data sets in order 
to develop studies of inter-annual variability [17]. 

Fig. 1 displays a Hovmöller plot [21] with the time evolution of 
CM SAF daily irradiation for 2010 and 2011 averaged along 10°W 
to 5°E, from 35.5°N to 44°N. Fig. 2 displays the average of annual 
global irradiation on the horizontal plane for 2010 and 2011. 


2.2. SIAR 


Land-measured daily irradiation is collected from the Agrocli- 
matic Information System for Irrigation (SIAR) [24] a free- 
download database operating since 1999, covering the majority 
of the irrigated area of Spain [25-29]. This network belongs to the 
Ministry of Agriculture, Food and Environment of Spain, as a tool 
to predict and study meteorological variables for agriculture. SIAR 
is composed by 12 regional centers and a national center, aiming 
to centralise and depurate measurements from the 361 stations of 
the network. The stations include SKYE-SP1110 (Campbell-Scien- 
tific)? or CMP6 (KIPP&ZONEN),? first class pyranometers according 
to the World Meteorological Organization (WMO). The absolute 
accuracy is within +5% and is typically lower than + 3%. 


2 ftp://ftp.campbellsci.com/pub/csl/outgoing/uk/manuals/sp1110.pdf. 
3 http://www.kippzonen.com/?product/1251/CMP-+ 6.aspx. 
4 http://www.wmo.int/pages/index_es.html. 


The calibration of the pyranometers is performed by Tragsatec 
[24,30] according to ISO 9847:1992 [31] using two CMP6- 
KIPP&ZONEN reference pyranometers [32] on a yearly basis. 
Irradiation is computed on a half-hourly basis from irradiance 
samples recorded each 10s, collated through a CR10X (Campbell 
Scientific) datalogger within the station and then sent to the 
regional and national centers [24]. 

Data has been filtered under two assumptions: average annual 
irradiation must be higher than 1000 kWh/m7, and only stations 
with more than 600 measurement days available (out of a total of 
730, 2 years) are selected. Besides, some stations have been 
omitted due to difficulties in the access to the coordinates of 
some stations, to uncompleted or spurious data series, or to 
stations out of the area of study. Eventually, 301 meteorological 
stations? (Fig. 3) and their daily global irradiation measurements 
on the horizontal plane for 2010 and 2011 have been considered. 


3. Methods 
3.1. Statistics 


The analysis is built upon the next structure (Fig. 4): 


e Analysis and comparison of daily and yearly global irradiation 
on the horizontal plane data obtained by on-ground 


5 The name and location data of these stations are available at http://solar. 
r-forge.r-project.org/data/SIAR.csv. 
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measurements, G*48(0), and satellite data, G°“"(0). The dif- 
ference between these sources is a matrix of values examined 
in Section 4.1. 


AGO) = G™S4F (Q)—GSAR(0) (1) 


e Analysis and comparison of yearly global irradiation incident 
on different planes (fixed, two axis, north-south horizontal 
axis) estimated from these two data sources, G3/"* and GgM*", 
respectively. The difference between these results is a set of 
three matrices examined in Section 4.1 


AGer = Gr = Gof (2) 


e Elaboration of a new map of yearly values with a smooth 
combination of both sources using geostatistical techniques, 
both for the horizontal plane, G? (0), and for the inclined 
planes, Gé?. Section 3.3 outlines the geostatistical interpola- 
tion technique (kriging with external drift, KED) used to 
combine the information from SIAR and CM SAF. 

The difference between the SIAR stations and the results of the 


interpolation are 


AG yep (0) = G*?(0)—-G>*"**(0) (3) 
AGF = GFF? — Cy" 4 


The difference between the CM SAF maps and the results of the 
interpolation are 


AGD (0) = G*P(0)—G™F(0) © 
ACH _ GKED _ GMSAF © 


These differences are summarized using several statistics: the 
unbiased and biased Root Mean Square Difference (RMSD* and RMSD, 
respectively), the Mean Bias Difference (MBD) and the Mean Absolute 
Difference (MAD) (Tables 3-5). It must be noted that RMSD? = 
RMSD*? + MBD? and that MAD < RMSD < n!/? . MAD. The reader is 
referred to Ref. [33] for the details on the convenience of these 
statistics. 

These statistics are normalized by the average SIAR values 
when comparing CM SAF or KED with SIAR (Eqs. (1)-(4)) and by 
the average CM SAF values when comparing KED with CM SAF 
(Eqs. (5) and (6)). For example, the RMSD*, RMSD, MBD and MAD 
corresponding to Eq. (1) are 


( [(c™ 0-670) = (o-c 0))] ‘ 1/2 
RMSD% = 


GAR (0) 
(7) 


Table 1 


(AG?(0)) we 
RMSDco = ~————— (8) 
G'AR (0) 
GAR (0) 
AG(0) 
MBDco = Zito, (10) 


3.2. Effective irradiation 
Three different tracking methods have been considered 


e Fixed plane, oriented towards the south and with optimum 
inclination angle (latitude minus 10°). 

e North-south horizontal axis tracker: the axis of rotation is 
horizontal with respect to the ground and is on a north-south 
line. Panels are mounted horizontally upon the tube which will 
rotate on its axis to track the apparent motion of the sun 
through the day. 

e Two-axis tracker: both azimuth and altitude are constantly 
changing to track the sun. 


Detailed description of these methods can be found in [34]. 
Table 1 summarizes the calculation procedure from global daily 
irradiation on the horizontal plane to effective global irradiation 
incident on an inclined plane. It is important to highlight that this 
calculation procedure does not include shadow losses. 

The first step of the procedure (once sun and trackers geome- 
try equations have been computed) is to decompose the daily 
global irradiation on the horizontal plane in two components, 
direct and diffuse irradiation. Diffuse fraction, the ratio of diffuse 
to global irradiation, is estimated from the clearness index with 
the equations proposed in [35]. 

The second step is to build irradiance profiles from daily 
irradiation values. The ratio of the diffuse irradiance to diffuse 
irradiation is assumed to be equivalent to the ratio of extra- 
terrestrial irradiance to extraterrestial irradiation. The ratio of 
global irradiance to daily global irradiation is estimated with 
the equations proposed in [35]. It must be noted that, because of 
the frequent low variability of solar irradiance, this step assumes 
that the average value of irradiance during a short time interval 
(for example, an hour) coincides numerically with the irradiation 
during that interval. Under this assumption the profile of irra- 
diance incident on a surface estimated in the next step can be 
aggregated to produce daily irradiation. 

The third step computes direct and diffuse irradiance incident 
on the inclined plane considering purely geometrical criteria. 
Direct irradiance is estimated with the solar zenith angle and 
the angle of incidence on the generator. Diffuse irradiance is 
calculated with the anisotropic model proposed in [36]. This 


Calculation procedure for the estimation of effective irradiation incident on a PV generator from daily global horizontal irradiation data. 


Step Method 


Sun and trackers geometry 

Decomposition of daily global horizontal 
irradiation 

Estimation of irradiance 


Set of equations as provided in [34] 


Correlation between diffuse fraction of horizontal irradiation and clearness index [35] 


Ratio of global irradiance to daily global irradiation [35] 


Estimation of irradiance on inclined surface The direct irradiance is calculated with geometrical equations. The estimation of the diffuse component makes use of the 


anisotropic model [36] 
Albedo irradiance 
Effects of dirt and angle of incidence 


Isotropic diffuse irradiance with reflection factor equal to 0.2. 
Equations proposed in [37]. A low constant dirtiness degree has been supposed (2%) 
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model divides the diffuse irradiance in circumsolar (near the sun 
region) and isotropic, using an anisotropy index to estimate the 
ratio between them. The albedo is assumed to be isotropic and is 
estimated from the global irradiance with a reflection factor 
of 0.2. 

The last step estimates the effective irradiance incident on a 
generator subtracting dust and angle of incidence losses from the 
incident irradiance with the model proposed in [37]. 


3.3. Geostatistical interpolation 


Geostatistics deals with the analysis of random fields Z(s), 
where s = (x,y) is a location and x, y its geographical coordinates. 
Measurements of the random field Z is commonly only available 
at a limited set of locations (in our case, the meteorological 
stations of the SIAR network). In order to predict the value of Z at 
locations without observations, the geostatistical analysis 
involves estimation and modelling of spatial correlation under 
simplifying assumptions of stationarity [38]. 

Assuming that the samples are representative, non- 
preferential and consistent, values of the field at a new location 
can be derived using a spatial prediction model. This geostatistical 
interpolation procedure is generally known as kriging [39]. A 
standard version of kriging is called ordinary kriging (OK). Here 
the predictions are based on the model: 


Z(S) = U+€(S) (11) 


where u is the constant stationary function (global mean) and €(s) 
is the spatially correlated stochastic part of variation. 

The assumption of constant mean is hardly acceptable for 
the estimation of irradiation over a large area. Ordinary kriging 
was initially tried within this study, generating inaccurate esti- 
mations due to the long distances among some stations. In 
mountainous-heterogeneous regions such as Galicia (north of 
Spain), this inaccuracy was more significant than in flat- 
homogeneous regions, such as Castilla-La-Mancha (center of 
Spain). 

This model can be improved including additional information 
from an exhaustively-sampled explanatory variable. If the expla- 
natory variable is significantly correlated with the field Z(s), 
predictions at a new location, sg, can be obtained modelling the 
deterministic and stochastic components separately 


Z(Sg) = M(Sg) + €(Sp) (12) 


where (Sg) is the value of the fitted deterministic part at the 
new location, €(Sg) is the interpolated residual. These two com- 
ponents can be derived with 


i n 
280) = D> BkaeSo)+ X Ai€(si) (13) 
k=0 i= 
where f, are the estimated coefficients of the deterministic 
model, q,(S) are the auxiliary predictors obtained from the fitted 
values of the explanatory variable at the new location, 4; are the 
kriging weights determined by the spatial dependence structure 
of the residual, and e(s;) are the residual at location s;. 

This improved model (Eq. (13)) is known as kriging with 
external drift (KED) or regression kriging [39]. In this paper, the 
explanatory variable is the irradiation on the horizontal plane 
estimated by CM SAF, G™^(0), both for the irradiation on the 
horizontal plane and for the irradiation incident on inclined 
planes. Therefore, the KED method is fed with three sources of 
information to produce new maps: 


e Yearly irradiation measurements on the horizontal plane from 
SIAR stations or estimations of yearly irradiation on the 
inclined plane based on the measurements from SIAR. 


Table 2 
Parameters of the variograms fitted to the SIAR data for different planes using CM 
SAF irradiation as explanatory variable. 


Irradiation plane Model Nugget Sill Range 
GO) Pure nugget 4609.11 - - 
Fixed Pure nugget 7275.30 - - 
N-S horizontal Spherical 13,138.08 6768.78 458.45 
Two axis Spherical 19,831.10 9336.59 478.47 
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Fig. 5. Semivariances and variograms fitted to the SIAR data for different planes 


using CM SAF irradiation as explanatory variable. 


e Estimations of yearly irradiation on the horizontal plane from 
CM SAF as explanatory variable. 

e A semivariogram function to model the spatial dependence 
structure of the residuals. 


The semivariogram is a function defined as [40,41] 
1 
7(h) = 5E(e(s)—e(s+h))° (14) 


where h is the separation vector between two locations, h = s;—s;. 
This equation is defined under the assumption that the variance 
of € is constant and that spatial correlation of € does not depend 
on location s but only on separation distance h. The estimator of 
the variogram, called the sample semivariogram is 


x 1 
ĵh) = IN, 280-8)" (15) 


with Np = {(S;,S;) : s—s;= h}, the set of all pairs of locations 
separated by vector h. It is common to assume that the variogram 
is isotropic and, consequently, that the correlation at two loca- 
tions depends only on the distance between them and not on the 
direction between them. 

The sample variogram gives estimates only at observed spatial 
lags. Therefore, it is not enough for prediction at new locations. A 
common solution is to infer a parametric variogram model from 
the data fitting a model to the sample variogram. Some well- 
known parametric variogram functions are the exponential, 
gaussian or spherical models. The parameters of the model to 
be determined are the sill, the range and the nugget [39]. Table 2 
displays the parameters of the variograms fitted to the SIAR data 
for different planes using CM SAF irradiation as explanatory 
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Fig. 6. Relative differences of the yearly values of horizontal (Eq. (1)) and effective irradiation (Eq. (2)) between CM SAF and SIAR for the whole set of SIAR stations. 
(a) Dotplot and (b) Map. 


Table 3 
Statistics of the yearly irradiation values from CM SAF and SIAR. The RMSD, RMSD*, MBD and MAD statistics are calculated with adimensionalized differences using GSO) 
or GJR as normalization factors. 
Irradiation plane Gcmsar (KWh/m?) Gsiar (KWh/m?) MBD (%) RMSD* (%) RMSD (%) MAD (%) 
GO 92.50 102.09 3.41 4.44 5.60 4.19 
Fixed 88.16 112.09 3.59 5.21 6.33 4.69 
N-S horiz 146.21 170.28 4.24 5.93 7.30 5.57 


Two-axis 155.99 195.63 4.33 6.36 7.69 5.84 
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Fig. 7. Time evolution of the relative differences between the daily global 
irradiation on the horizontal plane from SIAR and CM SAF. The red line represents 
the median and the blue lines represent the 5% and 95% quantiles. 
(For interpretation of the references to color in this figure caption, the reader is 
referred to the web version of this article.) 
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variable. Fig. 5 shows the semivariances and the fitted variogram 
models. 

The nugget effect, associated to micro-variability and mea- 
surement error, models the discontinuity of the variogram at the 
origin. When the nugget effect is present, the kriging method is 
not an exact interpolator (it does not preserve the original 
observations). It must be highlighted that the variograms corre- 
sponding to irradiation on the horizontal plane and on a fixed 
plane are the pure nugget model, that is, the residuals show no 
spatial auto-correlation. 


4. Discussion of the results 


4.1. Comparison between SIAR and CM SAF 


The comparison of G™^F (0) and G% (0) must be performed 
taking into account that the SIAR pyranometers present a toler- 
ance of 5% (Section 2.2). In Fig. 6a 71% of the locations are inside 
the range of this pyranometer uncertainty. Outside this 5% band, 
96.5% of the stations SIAR provide lower global irradiation values 
than CM SAF. 

The relative difference increases when a tracking system is 
considered, |AG(0)|/Gsiar(0) < |AGer|/Ger sian (Fig. 6). Besides, 
Table 3 shows that both the standard deviation of the irradiation 
values (@sjaz and Ocysar) and the statistics of the differences 
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Fig. 8. Statistics of the daily global irradiation on the horizontal plane from SIAR and CM SAF. (a) Absolute value of the MBD. (b) MAD. (c) RMSD and (d) RMSDc. 
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Fig. 9. Global solar irradiation estimated with KED using CM SAF as external drift. (a) GO. (b) Fixed. (c) NS Horiz. (d) Two. 
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Fig. 10. Histograms of the normalized differences between the effective irradia- 
tion incident on a fixed plane, a north-south horizontal axis tracker, and a two- 
axis tracker. 


Table 4 
Statistics of the yearly horizontal irradiation values from KED and SIAR. The RMSD, 
RMSD*, MBD and MAD statistics are calculated with adimensionalized differences 


using G**(0) or GF as normalization factors. 


efy 
Irradiation OKED OsIAR MBD  RMSD* RMSD MAD 
plane (kWh/m?) (kWh/m?) (%) (%) (%) (%) 
GO 72.67 102.09 0.00 4.28 4.28 3.02 
Fixed 63.01 112.09 0.00 5.23 5.23 3.66 
N-S horiz 122.67 170.28 0.00 4.60 4.60 3.32 
Two-axis 131.11 195.63 0.00 5.01 5.01 3.63 
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Fig. 11. Relative differences (%) between the horizontal irradiation estimated with 
kriging of the values at the SIAR stations using the Gy(0) from CM SAF as external 
drift, and the horizontal irradiation estimated from CM SAF. Positive values mean 
that the estimation with kriging is higher than with CM SAF. 


(RMSD and MAD) increase with the application of the formulas to 
account for tilted surfaces. This observation is consistent with 
[42,43], where the variability of the effective irradiation incident 
on tracking planes was reported to be higher than the variability 
of irradiation on the horizontal plane. 

No significant latitudinal behavior is appreciated in any of the 
cases of Fig. 6, although as per Figs. 1 and 2, solar irradiation is 
clearly latitudinally dependent. 

In Fig. 7, AG(O) presents a seasonal periodicity of the 5% and 95% 
quantiles, with a wider range for winter and more confined in 
summer. In this figure SIAR presents a set of samples in which 
GAR(Q) is significantly lower than G™^ (0). It may be explained 
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Fig. 12. Relative differences of the yearly values of horizontal (Eq. (3)) and effective irradiation (Eq. (4)) between KED and SIAR for the whole set of SIAR stations. 
(a) Dotplot. (b) Map. 
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Fig. 13. Relative differences (%) between the effective irradiation incident on a 
fixed plane estimated with kriging of the values at the SIAR stations using the 
G,(0) from CM SAF as external drift, and the effective irradiation estimated with 
the CM SAF raster. Positive values mean that the estimation with kriging is higher 
than with CM SAF. 
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Fig. 14. Relative differences (%) between the effective irradiation incident on a 
north-south horizontal axis tracker estimated with kriging of the values at the 
SIAR stations using the G,(0) from CM SAF as external drift, and the effective 
irradiation estimated with the CM SAF raster. Positive values mean that the 
estimation with kriging is higher than with CM SAF. 


due to local events not registered by the satellite resolution, or to 
failures in the on-ground registers, which were not detected when 
filtering spurious data. It is important to highlight that these 
extreme events are smoothed with the averages of annual sums. 
In Fig. 8a and b, the statistics MBD and MAD are lower than 5% 
in most of the stations, although a set of outliers is appreciated in 
the Valencia region (middle east of Spain). In Fig. 8c and d, the 
RMSD and RMSD* are generally lower than 7%. In a set of stations in 
the north of Spain in which the MBD were lower than 6%, the RMSD 
and RMSD* are significantly higher, which may be explained due to 
the strong meteorological variability existing in the Ebro valley. 
In the middle Ebro valley there are marked thermal contrasts, 
with possible generation of orographic fog, typical in valleys. 
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Fig. 15. Relative differences (%) between the effective irradiation incident on a 
two-axis tracker estimated with kriging of the values at the SIAR stations using the 
G,(0) from CM SAF as external drift, and the effective irradiation estimated with 
the CM SAF raster. Positive values mean that the estimation with kriging is higher 
than with CM SAF. 


This variability generates a much more variable distribution of 
error magnitudes which produces higher levels of RMSD [33]. 

Both Figs. 6 and 8 compare irradiation on the horizontal plane 
from SIAR and CM SAF. However, there are remarkable differences 
between them. For example, there are some stations in the north 
of Spain clearly visible in Fig. 6 (important difference between CM 
SAF and SIAR) but they are invisible in Fig. 8. To explain this 
apparently contradictory behavior is important to note that some 
stations include missing values in their data sets. Fig. 8 compares 
daily values with a collection of statistics computed without 
those missing values. However, Fig. 6 compares yearly values 
with missing values contributing as zeros. Therefore, those sta- 
tions with a higher proportion of missing values will provide 
lower annual irradiation values, although their daily statistics 
could be assimilable to a station without missing values. 


4.2. Comparison between KED, SIAR and CM SAF 


The KED technique does not perform as an interpolation 
function when the nugget effect is not null, which occurs in this 
study as shown in Table 2. This fact indicates that there is an 
intrinsic variability independent from the distance between sta- 
tions. In this case, the KED behaves as a smoothing function of the 
SIAR values, with the external drift of CM SAF, generating a 
solution that differs both from SIAR and CM SAF in the positions 
of the meteorological stations. 

Maps of global irradiation obtained with KED using CM SAF as 
external drift are shown in Fig. 9 on the horizontal surface, fixed 
tilted plane, tracking system with north-south axis and two-axis 
tracking system. Differences between one axis tracking and fixed 
plane range from 0.2% to 36%, and between two axis tracking and 
fixed plane range from 11% to 55% (Fig. 10). These differences are 
more significant in southern Huesca (42°N, 0°W), Zamora (42°N, 
6°W), the peninsular center (38°N to 41°N, 1°W to 6°W) and 
Almeria (37°N, 2°W), and lower along the Cantabric coast (43°N, 
1°W to 9°W), due to the reduced influence of direct irradiation. 
Differences in irradiation estimated with KED between the two 
tracking systems range from 11% to 14%, with higher values in the 
Ebro valley (40°N to 42°N, 5°W to 2°E) and along the Mediterra- 
nean coast and lower values in Jaen (38°N, 4°W). 
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Fig. 16. Violin plot of the relative differences (%) between irradation estimated with kriging of the values at the SIAR stations using CM SAF as external drift, and the 
irradiation estimated with the CM SAF raster. Positive values mean that the estimation with kriging is higher than with CM SAF. Each latitude interval include between 41 


and 43 stations with an overlap of 10%. 


Table 5 
Statistics of the yearly horizontal irradiation values from KED and CM SAF (Figs. 11 
and 13-15). The RMSD, RMSD*, MBD and MAD statistics are calculated with 


adimensionalized differences using G,“"(0) or Gg?" as normalization factors. 
OKED OCMSAF MBD RMSD* RMSD MAD 

(kWh/m?) (kWh/m?) (%) (%) (%) (%) 
GO 112.24 142.86 — 2.52 1.84 3.12 2.86 
Fixed 93.06 148.43 — 2.31 2.93 3.73 3.42 
N-S horiz 215.48 239.73 —3.48 2.62 4.35 3.60 
Two-axis 242.62 265.95 — 3.67 2.75 4.59 3.82 


Figs. 11, 6b and 12b reveal that |AGA(0)|/Gsiar(0) and 
[AGF ken | /Gep,siar are slightly lower than |AG(0)|/Gsar(0) and 
|AGer|/Ger,siar, respectively, and correspondingly for the CM SAF 
values. Once again, the RMSD and MAD values are very similar for 
all trackers and higher than those corresponding to the irradiation 
on the horizontal plane (Table 4). 

In Fig. 16,° higher latitudes present higher dispersion of the 
differences than lower latitudes, although values remain in a 4% 
band. Specially, from 40°N to north, just when average elevation 
increases (Fig. 3), dispersion values are higher. As already men- 
tioned, CM SAF shows a more inaccurate behavior when clouds or 
snow can appear. This fact can widen the range of differences for 
mountainous areas. In Fig. 2 mountainous areas act as modulators 
of irradiation [14]. In Fig. 11 (irradiation on the horizontal plane) 
SIAR only presents higher irradiation than CM SAF in the very 
north of Spain. Nevertheless, the variability of the previous map is 
in the range of 5%, which stands within the uncertainty band. 


6 This figure displays the data distribution with a violin plot, a combination of 
a boxplot and a kernel density plot. Therefore, this graphical tool shows the lower 
quartile, median (Q2), and the upper quartile, and the kernel density estimation. 


In Fig. 16, relative differences of irradiation incident on tilted 
planes reach values of 10% for fixed systems and — 10% for one- 
axis and two-axis with higher dispersion in these last cases. KED 
shows higher values than CM SAF especially in the north area for 
fixed systems, and to a lesser extent, for one and two axis around 
the Pyrenees area. The RMSD and MAD values are very similar for 
all trackers and higher than those corresponding to the irradiation 
on the horizontal plane (Table 5). 

One possible explanation for the positive values of relative 
differences existing in the area of Pyrenees would come from the 
influence of the terrain elevation on satellite methods [44]. The 
solar irradiation dependence with altitude is not well described in 
the satellite retrieving methods yet. In a mountain area each pixel 
cover an area of very varying altitude and therefore the irradia- 
tion estimations have more uncertainty than in flat terrains. 
Besides, satellite methods have difficulties in distinguishing snow 
and ice from clouds typically at mountain areas. An underestima- 
tion of the amount of water vapor and/or aerosols produces an 
overestimation in the values of CM SAF. Also the uncertainties 
associated with estimating the clouds cover plays an important 
role in divergence observed [14]. A special version of CM SAF for 
mountainous areas is being developed using methods proposed 
by and maybe, after this update, relative differences at those areas 
will decrease [45]. 


5. Conclusions 


An analysis and comparison of daily and yearly solar irradia- 
tion from the satellite CM SAF database and a set of 301 stations 
from the Spanish SIAR network is performed using data of 2010 
and 2011. This analysis is completed with the comparison of the 
estimations of effective irradiation incident on three different 
tilted planes (fixed, two axis tracking, north-south horizontal 
axis) using irradiation from these two data sources. Finally, a 
new map of yearly values of irradiation both on the horizontal 
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plane and on inclined planes is produced mixing both sources 
with geostatistical techniques (kriging with external drift, KED). 

The comparison between the irradiation values from SIAR, CM 
SAF and KED is performed in the context of the SIAR pyran- 
ometers tolerance (5%). The difference of global irradiation from 
SIAR and CM SAF at 71% of the locations are inside the range of 
this pyranometer uncertainty. Outside this 5% band, 96.5% of the 
SIAR stations provide lower global irradiation values than CM SAF. 
The relative difference increases when a tracking system is 
considered: both the standard deviation of the irradiation values 
(Osiar and Ocysar) and the statistics of the differences (RMSD and 
MAD) increase with the use of tracking systems. 

The Mean Absolute Difference (MAD) between CM SAF and 
SIAR is approximately 4% for the irradiation on the horizontal 
plane and is comprised between 5% and 6% for the irradiation 
incident on the inclined planes. 

The use of kriging with external drift reduces the difference 
with SIAR and CM SAF, both for the horizontal plane and for 
inclined planes. The MAD between KED and SIAR, and KED and 
CM SAF is approximately 3% for the irradiation on the horizontal 
plane and is comprised between 3% and 4% for the irradiation 
incident on the inclined planes. 


6. Software 


The methods described in this paper have been implemented 
using the free software environment R [46] and several contrib- 
uted packages, namely: gstat [40] and sp [41] for the geosta- 
tistical analysis; solar [47] for the solar geometry, irradiation 
and PV energy calculations; raster [48] for spatial data manip- 
ulation and analysis, and rasterVis [49] for spatial data visua- 
lization methods. 

The source code is available at https://github.com/oscarperpi 
nan/CMSAF-SIAR. 
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